Plants.f90 Source File
Simulate vegetation dynamic
Source Code
!! Simulate vegetation dynamic !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.0 - 31st October 2018 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 31/Oct/2018 | Original code | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! Module to model vegetation dynamic ! MODULE Plants ! ! Modules used: USE DataTypeSizes, ONLY : & ! Imported Type Definitions: short, float USE LogLib, ONLY: & !Imported routines: Catch USE GridLib, ONLY: & !imported definitions: grid_real, grid_integer, & !imported routines: NewGrid, GridDestroy, & !imported parameters: NET_CDF USE GridOperations, ONLY : & !imported routines: GridByIni, CellArea USE Chronos, ONLY: & !imported definitions: DateTime, & !Imported routines: DayOfYear, GetYear, & !Imported operations: OPERATOR (==), & OPERATOR ( > ) USE IniLib, ONLY : & !Imported derived types: IniList, & !Imported routines: IniOpen, IniClose , & IniReadInt , IniReadReal, & IniReadString, SectionIsPresent, & KeyIsPresent USE StringManipulation, ONLY : & !imported routines ToString USE PlantsAllometrics, ONLY : & !imported routines: CrownDiameter, CanopyCover USE PlantsModifiers, ONLY: & !imported routines: AGEmod, TEMPmod, SWCmod, VPDmod, CO2mod USE PlantsInterception, ONLY: & !imported variables interceptionParametersByMap, & laimax, canopymax USE PlantsManagement, ONLY: & !imported definition: Practice, & !imported variables: plants_management, & management_map, & !imported routines: SetPlantsManagement, & ApplyPlantsManagement, & SetPractice USE PlantsMortality, ONLY : & !imported routines: KillPlants USE Units, ONLY : & ! imported parameters: hectare, pi IMPLICIT NONE ! Global (i.e. public) Declarations: PUBLIC :: PlantsConfig PUBLIC :: PlantsGrow PUBLIC :: PlantsParameterUpdate INTEGER (KIND = short) :: dtPlants TYPE (grid_real) :: lai !!leaf area index (m2/m2) TYPE (grid_real) :: fvcover !!fraction of cell covered by vegetation TYPE (grid_real) :: rsMin !!minimum stomatal resistance TYPE (grid_real) :: gpp !! carbon gross primary production (t) TYPE (grid_real) :: npp !! carbon net primary production (t) TYPE (grid_real) :: carbonroot !! carbon mass of root (t) TYPE (grid_real) :: carbonstem !! carbon mass of stem (t) TYPE (grid_real) :: carbonleaf !! carbon mass of leaf (t) TYPE (grid_real) :: dbh !!diameter at brest height (cm) TYPE (grid_real) :: plantsHeight !! tree height (m) TYPE (grid_real) :: density !! Trees/hectare TYPE (grid_real) :: stemyield !! stem yield (t) LOGICAL :: rsMinLoaded = .FALSE. LOGICAL :: fvcoverLoaded = .FALSE. LOGICAL :: laiLoaded = .FALSE. LOGICAL :: plantsHeightLoaded = .FALSE. LOGICAL :: updatePlantsParameters = .FALSE. TYPE PlantsSpecies !Descriptor CHARACTER (LEN = 50) :: name !!species name, Fagus, Abies, ... CHARACTER (LEN = 1) :: phenology !! evergreen (E), deciduous (D) REAL (KIND = float) :: alpha !! canopy quantum efficiency (molC/molPAR) REAL (KIND = float) :: GPPtoNPP !!GPP/NPP ratio REAL (KIND = float) :: k !! extinction coefficient for absorption of PAR by canopy REAL (KIND = float) :: cra !! Chapman-Richards asymptotic maximum height REAL (KIND = float) :: crb !! Chapman-Richards exponential decay parameter REAL (KIND = float) :: crc !! Chapman-Richards shape parameter REAL (KIND = float) :: as !! scaling coefficient in stem mass v diameter relationship REAL (KIND = float) :: ns !! scaling exponent in stem mass v diameter relationship !REAL (KIND = float) :: saf !! stem allocation factor !REAL (KIND = float) :: raf !! root allocation factor !REAL (KIND = float) :: laf !! leaf (foliage) allocation factor REAL (KIND = float) :: dbhdcmin !! minimum ratio between stem and crown diameters (m/cm) REAL (KIND = float) :: dbhdcmax !! maximum ratio between stem and crown diameters (m/cm) REAL (KIND = float) :: denmin !! minimum tree density (trees/ha) REAL (KIND = float) :: denmax !! minimum tree density (trees/ha) REAL (KIND = float) :: agemax !! maximum age REAL (KIND = float) :: tmax !! maximum temperature for vegetation growing (°C) REAL (KIND = float) :: tmin !! minimum temperature for vegetation growing (°C) REAL (KIND = float) :: topt !! optimum temperature for vegetation growing (°C) REAL (KIND = float) :: theta_fswc !! parameter to compute soil water content modifier REAL (KIND = float) :: theta_fvpd !! parameter to compute vapor pressure deficit modifier REAL (KIND = float) :: fpra !! parameter 1 to compute allocation factors REAL (KIND = float) :: fprn !! parameter 2 to compute allocation factors REAL (KIND = float) :: spra !! parameter 3 to compute allocation factors REAL (KIND = float) :: sprn !! parameter 4 to compute allocation factors REAL (KIND = float) :: ltr !!leaf turnover rate [s-1] REAL (KIND = float) :: rtr !!root turnover rate [s-1] REAL (KIND = float) :: tcold_leaf !! temperature threshold that accelerates leaf turnover (°C) REAL (KIND = float) :: sla !! specific leaf area [m2/Kg] REAL (KIND = float) :: hdmin !! H/D ratio in carbon partitioning for low density REAL (KIND = float) :: hdmax !! H/D ratio in carbon partitioning for high density REAL (KIND = float) :: albedo !!plant albedo (0-1) REAL (KIND = float) :: laimax !!maximum leaf area index used for precipitation interception (m2/m2) REAL (KIND = float) :: canopymax !! maximum canopy storage capacity (m) REAL (KIND = float) :: wood_density !! wood density (kg/m3) REAL (KIND = float) :: ms !! Fractions of mean stem biomass pools per tree on each dying tree REAL (KIND = float) :: mr !! Fractions of mean root biomass pools per tree on each dying tree REAL (KIND = float) :: mf !! Fractions of mean leaf biomass pools per tree on each dying tree REAL (KIND = float) :: wSx1000 !! Fractions of mean leaf biomass pools per tree on each dying tree END TYPE PlantsSpecies TYPE PlantsCohort !group of plants of the same species that are born during the same year TYPE (PlantsSpecies) :: species !state variables REAL (KIND = float) :: age !!age [years] REAL (KIND = float) :: gpp !! gross primary production (t/ha) REAL (KIND = float) :: npp !! net primary production (t/ha) REAL (KIND = float) :: lai !! leaf area index [m2/m2] REAL (KIND = float) :: mass_root !! mass of root (t/ha) REAL (KIND = float) :: mass_stem !! mass of stem (t/ha) REAL (KIND = float) :: mass_stem_year_previous !! mass of stem of the previous year (t/ha) REAL (KIND = float) :: mass_leaf !! mass of leaf (t/ha) REAL (KIND = float) :: mass_total !! total biomass root + stem + foliage (t/ha) REAL (KIND = float) :: stem_yield !! stem mass produced by plants cut (t/ha) REAL (KIND = float) :: density !! number of plants per hectar REAL (KIND = float) :: height !! plants heigth [m] REAL (KIND = float) :: dbh !!diameter at brest heigth [cm] REAL (KIND = float) :: canopy_cover !! percentage of surface stand covered by canopy [0-1] REAL (KIND = float) :: crown_diameter !! diameter of crown [m] REAL (KIND = float) :: apar !! absorbed photosynthetically active radiation (molPAR m-2) TYPE(PlantsCohort), POINTER :: next =>null() ! pointer to the next cohort in the same stand END TYPE PlantsCohort TYPE PlantsStand ! list of cohorts TYPE (PlantsCohort), POINTER, PRIVATE :: first => null() !TYPE (PlantsCohort), POINTER, PRIVATE :: last => null() !TYPE (PlantsCohort), POINTER, PRIVATE :: iter => null() INTEGER (KIND = short) , PRIVATE :: lenght = 0 INTEGER (KIND = short) , PRIVATE :: i !row of the stand in grid plants_mask INTEGER (KIND = short) , PRIVATE :: j !column of the stand in grid plants_mask TYPE (Practice) :: thinning END TYPE PlantsStand !TYPE PlantsAge ! REAL (KIND = float) :: value !!age [yesrs] ! INTEGER (KIND = short) :: period !! 0 for adult tree, 1 for young tree ! INTEGER (KIND = short) :: species_count !! number of species in the cell ! TYPE (PlantsSpecies), ALLOCATABLE :: species (:) !END TYPE PlantsAge ! !TYPE PlantsHeight ! REAL (KIND = float) :: value !!value of height class [m] ! REAL (KIND = float) :: layer_coverage !! cell coverage of that layer ! INTEGER (KIND = short) :: dominance !! -1 no trees in veg period, 1 trees in veg period ! INTEGER (KIND = short) :: ages_count !! counter of age classes in the cell ! TYPE (PlantsAge), ALLOCATABLE :: ages (:) !END TYPE PlantsHeight !local (i.e. private) declarations: !variables and parameters LOGICAL, PRIVATE :: simulatePlants != 1 when plants dynamic is simulated, = 0 read vegetation parameters from files LOGICAL, PRIVATE :: useCO2modifier ! define wheter use the CO2 modifier LOGICAL, PRIVATE :: mortality ! wheter to apply mortality or not TYPE (grid_integer), PRIVATE :: plants_mask !!define cells active for plants dynamic simulation TYPE (PlantsStand), PRIVATE, ALLOCATABLE :: forest (:) !vector of stands INTEGER (KIND = short), PRIVATE :: count_stands !number of stands in the forest (equal to the number of cells in the domain) TYPE (PlantsSpecies), ALLOCATABLE, PRIVATE :: species (:) ! !REAL (KIND = float), PRIVATE :: rs_par !! fraction of shortwave radiation to compute PAR [0-1] INTEGER (KIND = short), PRIVATE :: year_prev, year_new REAL (KIND = float), PRIVATE, PARAMETER :: C_molar_mass = 0.0120107 ! kg/mol !routines PRIVATE :: AparCalc PRIVATE :: ReadSpecies PRIVATE :: SetIC PRIVATE :: GrowLeaf PRIVATE :: GrowDBH PRIVATE :: ForestToGrid !======= CONTAINS !======= !============================================================================== !| Description: ! Initialize Plants module SUBROUTINE PlantsConfig & ! (iniFile, mask, begin, end) IMPLICIT NONE ! Arguments with intent(in): CHARACTER (LEN = 300), INTENT(IN) :: iniFile !!configuration file TYPE (grid_integer), INTENT(IN) :: mask !!mask of simulation domain TYPE (DateTime), INTENT (IN) :: begin !!simulation starting date TYPE (DateTime), INTENT(IN) :: end !!simulation ending date !local declarations: TYPE(IniList) :: iniDB !!store configuration info INTEGER (KIND = short) :: i, j, k CHARACTER (LEN = 300) :: species_file CHARACTER (LEN = 300) :: ICfile CHARACTER (LEN = 300) :: PMfile REAL (KIND = float) :: scalar !------------end of declaration------------------------------------------------ CALL IniOpen (iniFile, iniDB) !need plants dynamic simulation? simulatePlants = IniReadInt ('plants-simulation', iniDB) IF ( simulatePlants == 1) THEN !plants dynamic is simulated: configure forest model !plants mask IF ( SectionIsPresent ( 'plants-mask', iniDB) ) THEN CALL GridByIni (iniDB, plants_mask, section = 'plants-mask') ELSE CALL Catch ('error', 'Plants', 'Mask missing in configuration file') END IF !CO2 modifier IF ( SectionIsPresent ( 'use-co2-modifier', iniDB) ) THEN IF ( IniReadInt ('use-co2-modifier', iniDB) == 1 ) THEN useCO2modifier = .TRUE. ELSE useCO2modifier = .FALSE. END IF ELSE !if key is not found, set it to default false useCO2modifier = .FALSE. END IF !mortality IF ( KeyIsPresent ( 'mortality', iniDB) ) THEN IF ( IniReadInt ('mortality', iniDB) == 1 ) THEN mortality = .TRUE. ELSE mortality = .FALSE. END IF ELSE !if mortality is not found, assume it is not considered mortality = .FALSE. END IF !load species species_file = IniReadString ('plants-species-file', iniDB) CALL ReadSpecies (species_file) !compute number of stands and allocate count_stands = 0 DO i = 1, plants_mask % idim DO j = 1, plants_mask % jdim IF ( plants_mask % mat (i,j) /= plants_mask % nodata ) THEN count_stands = count_stands + 1 END IF END DO END DO ALLOCATE ( forest ( count_stands) ) !initialize list k = 1 DO i = 1, plants_mask % idim DO j = 1, plants_mask % jdim IF ( plants_mask % mat (i,j) /= plants_mask % nodata ) THEN !initialize the first cohort in the stand ALLOCATE ( forest (k) % first) NULLIFY ( forest (k) % first % next) forest (k) % lenght = 1 !store position of stand forest (k) % i = i forest (k) % j = j k = k + 1 END IF END DO END DO !set initial condition ICfile = IniReadString ('plants-ic-file', iniDB) CALL SetIC (ICfile) !set plants management options IF (KeyIsPresent ('plants-management-file',iniDB) ) THEN plants_management = .TRUE. PMfile = IniReadString ('plants-management-file', iniDB) CALL SetPlantsManagement (PMfile, begin, end) !set management pratice for each stand DO k = 1, count_stands CALL SetPractice ( management_map % mat (forest (k) % i, forest (k) % j) , forest (k) % thinning ) END DO ELSE !management is turned off plants_management = .FALSE. END IF !set interception maps if required IF (interceptionParametersByMap == 0 ) THEN CALL Catch ('warning', 'Plants', 'interception parameters are set only where plants map is not nodata') !allocate maps CALL NewGrid (laimax, plants_mask) CALL NewGrid (canopymax, plants_mask) DO k = 1, count_stands i = forest (k) % i j = forest (k) % j laimax % mat (i,j) = forest (k) % first % species % laimax canopymax % mat (i,j) = forest (k) % first % species % canopymax END DO END IF !allocate grids for state variables CALL NewGrid (lai, plants_mask) CALL NewGrid (fvcover, plants_mask) CALL NewGrid (gpp, plants_mask) CALL NewGrid (npp, plants_mask) CALL NewGrid (carbonroot, plants_mask) CALL NewGrid (carbonstem, plants_mask) CALL NewGrid (carbonleaf, plants_mask) CALL NewGrid (dbh, plants_mask) CALL NewGrid (plantsHeight, plants_mask) CALL NewGrid (density, plants_mask) CALL NewGrid (stemyield, plants_mask) !fill in grids CALL ForestToGrid (lai, 'lai') CALL ForestToGrid (fvcover, 'fv') CALL ForestToGrid (gpp, 'gpp') CALL ForestToGrid (npp, 'npp') CALL ForestToGrid (carbonroot, 'root') CALL ForestToGrid (carbonstem, 'stem') CALL ForestToGrid (carbonleaf, 'leaf') CALL ForestToGrid (dbh, 'dbh') CALL ForestToGrid (plantsHeight, 'height') CALL ForestToGrid (density, 'density') CALL ForestToGrid (stemyield, 'stemyield') fvcoverLoaded = .TRUE. laiLoaded = .TRUE. plantsHeightLoaded = .TRUE. ELSE ! no plants dynamic simulation required: read parameters from files !leaf area index IF ( SectionIsPresent ( 'lai', iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, 'lai') ) THEN scalar = IniReadReal ('scalar', iniDB, 'lai') CALL NewGrid (lai, mask, scalar) ELSE !check if parameter may change in time IF (KeyIsPresent ('format', iniDB, 'lai') ) THEN IF ( IniReadString ('format', iniDB, 'lai') == 'net-cdf' ) THEN updatePlantsParameters = .TRUE. END IF END IF !read map CALL GridByIni (iniDB, lai, section = 'lai') END IF laiLoaded = .TRUE. ELSE CALL Catch ('warning', 'Plants', 'LAI missing in configuration file') END IF !fraction of vegetation coverage IF ( SectionIsPresent ( 'vegetation-fraction', iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, 'vegetation-fraction') ) THEN scalar = IniReadReal ('scalar', iniDB, 'vegetation-fraction') CALL NewGrid (fvcover, mask, scalar) ELSE !check if parameter may change in time IF (KeyIsPresent ('format', iniDB, 'vegetation-fraction') ) THEN IF ( IniReadString ('format', iniDB, 'vegetation-fraction') == 'net-cdf' ) THEN updatePlantsParameters = .TRUE. END IF END IF !read map CALL GridByIni (iniDB, fvcover, section = 'vegetation-fraction') END IF fvcoverLoaded = .TRUE. ELSE CALL Catch ('warning', 'Plants', 'vegetation-fraction is missing') END IF !plants height IF ( SectionIsPresent ( 'vegetation-height', iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, 'vegetation-height') ) THEN scalar = IniReadReal ('scalar', iniDB, 'vegetation-height') CALL NewGrid (plantsHeight, mask, scalar) ELSE !check if parameter may change in time IF (KeyIsPresent ('format', iniDB, 'vegetation-height') ) THEN IF ( IniReadString ('format', iniDB, 'vegetation-height') == 'net-cdf' ) THEN updatePlantsParameters = .TRUE. END IF END IF !read map CALL GridByIni (iniDB, plantsHeight, section = 'vegetation-height') END IF plantsHeightLoaded = .TRUE. ELSE CALL Catch ('warning', 'Plants', 'vegetation-height is missing') END IF END IF !minimum stomatal resistance. this is not computed by forest model and it may be used for computing PET IF ( SectionIsPresent ( 'min-stomatal-resistance', iniDB) ) THEN IF (KeyIsPresent ('scalar', iniDB, 'min-stomatal-resistance') ) THEN scalar = IniReadReal ('scalar', iniDB, 'min-stomatal-resistance') CALL NewGrid (rsMin, mask, scalar) ELSE !check if parameter may change in time IF (KeyIsPresent ('format', iniDB, 'min-stomatal-resistance') ) THEN IF ( IniReadString ('format', iniDB, 'min-stomatal-resistance') == 'net-cdf' ) THEN updatePlantsParameters = .TRUE. END IF END IF !read map CALL GridByIni (iniDB, rsMin, section = 'min-stomatal-resistance') END IF rsMinLoaded = .TRUE. ELSE CALL Catch ('warning', 'Plants', 'min-stomatal-resistance is missing') END IF CALL IniClose (iniDB) RETURN END SUBROUTINE PlantsConfig !============================================================================== !| Description: ! update parameter map that change in time SUBROUTINE PlantsParameterUpdate & ! (time) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT (IN) :: time !local declarations: CHARACTER (LEN = 300) :: filename CHARACTER (LEN = 300) :: varname !------------------------------end of declarations----------------------------- !leaf area index IF ( time == lai % next_time ) THEN !destroy current grid filename = lai % file_name varname = lai % var_name CALL GridDestroy (lai) !read new grid in netcdf file CALL NewGrid (lai, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !fraction of vegetation coverage IF ( time == fvcover % next_time ) THEN !destroy current grid filename = fvcover % file_name varname = fvcover % var_name CALL GridDestroy (fvcover) !read new grid in netcdf file CALL NewGrid (fvcover, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !plants height IF ( time == plantsHeight % next_time ) THEN !destroy current grid filename = plantsHeight % file_name varname = plantsHeight % var_name CALL GridDestroy (plantsHeight) !read new grid in netcdf file CALL NewGrid (plantsHeight, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !minimum stomatal resistance IF ( time == rsMin % next_time ) THEN !destroy current grid filename = rsMin % file_name varname = rsMin % var_name CALL GridDestroy (rsMin) !read new grid in netcdf file CALL NewGrid (rsMin, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF RETURN END SUBROUTINE PlantsParameterUpdate !============================================================================== !| Description: ! Update plants state variables SUBROUTINE PlantsGrow & ! (time, radiation, temperature, swc, sfc, swp, rh, co2) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT(IN) :: time !!current time TYPE (grid_real), INTENT(IN) :: radiation !!shortwave radiation [w/m2] TYPE (grid_real), INTENT(IN) :: temperature !!air temperature [°C] TYPE (grid_real), INTENT(IN) :: swc !! soil water content [m3/m3] TYPE (grid_real), INTENT(IN) :: sfc !! soil field capacity [m3/m3] TYPE (grid_real), INTENT(IN) :: swp !! soil wilting point [m3/m3] TYPE (grid_real), INTENT(IN) :: rh !! air relative humidity [0-1] REAL (KIND = float), OPTIONAL, INTENT(IN) :: co2 !! CO2 [ppm] !local declarations: CHARACTER (LEN = 300) :: filename, varname INTEGER (KIND = short) :: k, l TYPE (PlantsCohort), POINTER :: cohort REAL (KIND = float) :: fSWC !! GPP soil water content modifier REAL (KIND = float) :: fAGE !! GPP age modifier REAL (KIND = float) :: fTEMP !! GPP temperature modifier REAL (KIND = float) :: fVPD !! GPP vapor pressure deficit modififer REAL (KIND = float) :: fCO2 !! GPP CO2 modifier !REAL (KIND = float) :: deltaMroot !! root mass increment !REAL (KIND = float) :: deltaMstem !! stem mass increment REAL (KIND = float) :: rootaf !!root allocation factor REAL (KIND = float) :: stemaf !! stem allocation factor REAL (KIND = float) :: leafaf !! leaf allocation factor REAL (KIND = float) :: pfs !!used to compute stem allocation factor LOGICAL :: is_new_year !! flag to mark new year has begun TYPE (Practice) :: practice !------------end of declaration------------------------------------------------ ! need to update vegetation parameters from file? IF ( simulatePlants == 0) THEN !check leaf area index IF ( time == lai % next_time ) THEN !destroy current grid filename = lai % file_name varname = lai % var_name CALL GridDestroy (lai) !read grid in netcdf file CALL NewGrid (lai, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !check vegetation fraction IF ( time == fvcover % next_time ) THEN !destroy current grid filename = fvcover % file_name varname = fvcover % var_name CALL GridDestroy (fvcover) !read grid in netcdf file CALL NewGrid (fvcover, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !check plants height IF ( time == plantsHeight % next_time ) THEN !destroy current grid filename = plantsHeight % file_name varname = plantsHeight % var_name CALL GridDestroy (plantsHeight) !read grid in netcdf file CALL NewGrid (plantsHeight, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !check minimum stomatal resistance IF ( rsMinLoaded .AND. time == rsMin % next_time ) THEN !destroy current grid filename = rsMin % file_name varname = rsMin % var_name CALL GridDestroy (rsMin) !read grid in netcdf file CALL NewGrid (rsMin, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF ELSE !simulate forest dynamics year_new = GetYear (time) DO k = 1, count_stands cohort => forest (k) % first DO l = 1, forest (k) % lenght !loop through all cohorts cohort % stem_yield = 0. !check management practices IF ( plants_management ) THEN IF ( ALLOCATED (forest (k) % thinning % cuts) ) THEN practice = forest (k) % thinning IF (time > practice % next .AND. practice % current < SIZE ( practice % cuts )) THEN practice % current = practice % current + 1 IF ( practice % cuts ( practice % current ) % reforestation ) THEN !set species parameters cohort % species = species ( practice % cuts ( practice % current ) % species) END IF CALL ApplyPlantsManagement (time, practice, cohort % density, & cohort % mass_root, cohort % mass_stem, & cohort % mass_leaf, cohort % mass_total, & cohort % lai, cohort % canopy_cover, & cohort % age, cohort % height, cohort % dbh, & cohort % stem_yield ) !update crown diameter cohort % crown_diameter = CrownDiameter ( dbh = cohort % dbh, & den = cohort % density, & denmin = cohort % species % denmin, & denmax = cohort % species % denmax, & dbhdcmin = cohort % species % dbhdcmin, & dbhdcmax = cohort % species % dbhdcmax ) !update canopy cover cohort % canopy_cover = CanopyCover ( cohort % crown_diameter, cohort % density ) !update time of next cut IF ( practice % current < SIZE ( practice % cuts ) ) THEN practice % next = practice % cuts ( practice % current + 1 ) % time END IF !copy practice back forest (k) % thinning = practice END IF END IF END IF !update age IF ( DayOfYear (time) == 1 .AND. year_new /= year_prev ) THEN is_new_year = .TRUE. cohort % age = cohort % age + 1 cohort % mass_stem_year_previous = cohort % mass_stem ELSE is_new_year = .FALSE. END IF !compute APAR cohort % apar = AparCalc (rad = radiation % mat (forest (k) % i, forest (k) % j) , & lai = cohort % lai, & k = cohort % species % k, & alb = cohort % species % albedo, & dt = dtPlants) !compute CO2 modifier IF ( useCO2modifier ) THEN fCO2 = CO2mod (co2, cohort % age ) ELSE fCO2 = 1. END IF !compute age modifier fAGE = AGEmod ( cohort % age, cohort % species % agemax) !compute temperature modifier fTEMP = TEMPmod (Ta = temperature % mat (forest (k) % i, forest (k) % j), & Tmin = cohort % species % Tmin, & Tmax = cohort % species % Tmax, & Topt = cohort % species % Topt ) !compute soil water content modifier fSWC = SWCmod (swc = swc % mat (forest (k) % i, forest (k) % j), & wp = swp % mat (forest (k) % i, forest (k) % j), & fc = sfc % mat (forest (k) % i, forest (k) % j), & theta = cohort % species % theta_fswc ) !compute vapor pressure deficit modifier fVPD = VPDmod (Ta = temperature % mat (forest (k) % i, forest (k) % j), & RH = rh % mat (forest (k) % i, forest (k) % j), & kd = cohort % species % theta_fvpd) !compute gross primary production cohort % gpp = cohort % apar * & ! absorbed photosynthetically active radiation molPAR m-2 cohort % species % alpha * &! canopy quantum efficiency (molC/molPAR) C_molar_mass * & !conversion to Kg hectare / & ! m2 in one hectare 1000. * & !conversion to t fAGE * & !age modifier fTEMP * & !temperature modifier fSWC * & !soil water content modifier fVPD * &!vapor pressure deficit modifier fCO2 !compute net primary production cohort % npp = cohort % gpp * cohort % species % GPPtoNPP !compute root allocation factor rootaf = 0.5 / ( 1. + 2.5 * fAGE * fTEMP * fSWC ) !compute stem allocation factor pfs = (cohort % species % fprn * cohort % species % fpra) / & (cohort % species % sprn * cohort % species % spra) * & (cohort % density * cohort % dbh / 100.) ** & (cohort % species % fprn - cohort % species % sprn) stemaf = ( 1. - rootaf ) / ( 1. + pfs ) !compute leaf allocation factor leafaf = 1. - rootaf - stemaf ! update stem biomass cohort % mass_stem = cohort % mass_stem + stemaf * cohort % npp !update root biomass cohort % mass_root = cohort % mass_root + & !current mass rootaf * cohort % npp - & !mass increment cohort % species % rtr * cohort % mass_root * dtPlants !turnover rate ! update leaf biomass and leaf area index CALL GrowLeaf (npp = cohort % npp , & af = leafaf, & Ta = temperature % mat (forest (k) % i, forest (k) % j), & Tcold = cohort % species % tcold_leaf, & swc = swc % mat (forest (k) % i, forest (k) % j), & swp = swp % mat (forest (k) % i, forest (k) % j), & sfc = sfc % mat (forest (k) % i, forest (k) % j), & tr = cohort % species % ltr, & sla = cohort % species % sla, & mleaf = cohort % mass_leaf, & lai = cohort % lai ) !update total biomass cohort % mass_total = cohort % mass_root + cohort % mass_stem + cohort % mass_leaf !update dbh and height tree every new year !IF (is_new_year) THEN !CALL GrowDBH (cc = cohort % canopy_cover, & ! hdmin = cohort % species % hdmin, & ! hdmax = cohort % species % hdmax, & ! ws = cohort % mass_stem, & ! dws = cohort % mass_stem - cohort % mass_stem_year_previous, & ! !dws = stemaf * cohort % npp, & ! DBH = cohort % dbh, & ! height = cohort % height) CALL GrowDBHech2o (cc = cohort % canopy_cover, & hdmin = cohort % species % hdmin, & hdmax = cohort % species % hdmax, & ws = cohort % mass_stem, & !dws = cohort % mass_stem - cohort % mass_stem_year_previous, & dws = stemaf * cohort % npp, & DBH = cohort % dbh, & height = cohort % height, & tree_density = cohort % density, & wood_density = cohort % species % wood_density, & age = cohort % age, & maxage = cohort % species % agemax ) !END IF !update crown diameter cohort % crown_diameter = CrownDiameter ( dbh = cohort % dbh, den = cohort % density, & denmin = cohort % species % denmin, & denmax = cohort % species % denmax, & dbhdcmin = cohort % species % dbhdcmin, & dbhdcmax = cohort % species % dbhdcmax ) !update canopy cover cohort % canopy_cover = CanopyCover ( cohort % crown_diameter, cohort % density ) !apply mortality IF (mortality) THEN CALL KillPlants (dtPlants, cohort % age, cohort % species % agemax, & cohort % species % ms, cohort % species % mf, & cohort % species % mr, cohort % species % wSx1000, & cohort % crown_diameter, cohort % density, & cohort % canopy_cover, cohort % mass_stem, & cohort % mass_root, cohort % mass_leaf, & cohort % mass_total) END IF !jump to next cohort cohort => cohort % next END DO END DO year_prev = year_new !update state variable grids CALL ForestToGrid (lai, 'lai') CALL ForestToGrid (fvcover, 'fv') CALL ForestToGrid (gpp, 'gpp') CALL ForestToGrid (npp, 'npp') CALL ForestToGrid (carbonroot, 'root') CALL ForestToGrid (carbonstem, 'stem') CALL ForestToGrid (carbonleaf, 'leaf') CALL ForestToGrid (dbh, 'dbh') CALL ForestToGrid (plantsHeight, 'height') CALL ForestToGrid (density, 'density') CALL ForestToGrid (stemyield, 'stemyield') END IF RETURN END SUBROUTINE PlantsGrow !============================================================================== !| Description: ! read species parameters SUBROUTINE ReadSpecies & ! (inifile) IMPLICIT NONE CHARACTER (LEN = *), INTENT(in) :: inifile !local declarations: TYPE (IniList) :: speciesDB INTEGER (KIND = short) :: nspecies INTEGER (KIND = short) :: i !------------------------------end of declarations---------------------------- ! open and read configuration file CALL IniOpen (inifile, speciesDB) nspecies = IniReadInt ('species', speciesDB) ALLOCATE (species (nspecies) ) DO i = 1, nspecies species (i) % name = IniReadString ('name', speciesDB, section = ToString(i)) species (i) % phenology = IniReadString ('phenology', speciesDB, section = ToString(i)) species (i) % alpha = IniReadReal ('alpha', speciesDB, section = ToString(i)) species (i) % GPPtoNPP = IniReadReal ('GPP-NPP', speciesDB, section = ToString(i)) species (i) % k = IniReadReal ('k', speciesDB, section = ToString(i)) species (i) % cra = IniReadReal ('cra', speciesDB, section = ToString(i)) species (i) % crb = IniReadReal ('crb', speciesDB, section = ToString(i)) species (i) % crc = IniReadReal ('crc', speciesDB, section = ToString(i)) species (i) % as = IniReadReal ('as', speciesDB, section = ToString(i)) species (i) % ns = IniReadReal ('ns', speciesDB, section = ToString(i)) !species (i) % saf = IniReadReal ('stem-allocation-factor', speciesDB, section = ToString(i)) !species (i) % raf = IniReadReal ('root-allocation-factor', speciesDB, section = ToString(i)) !species (i) % laf = IniReadReal ('leaf-allocation-factor', speciesDB, section = ToString(i)) species (i) % dbhdcmin = IniReadReal ('dbhdcmin', speciesDB, section = ToString(i)) species (i) % dbhdcmax = IniReadReal ('dbhdcmax', speciesDB, section = ToString(i)) species (i) % denmin = IniReadReal ('denmin', speciesDB, section = ToString(i)) species (i) % denmax = IniReadReal ('denmax', speciesDB, section = ToString(i)) species (i) % agemax = IniReadReal ('agemax', speciesDB, section = ToString(i)) species (i) % tmin = IniReadReal ('tmin', speciesDB, section = ToString(i)) species (i) % tmax = IniReadReal ('tmax', speciesDB, section = ToString(i)) species (i) % topt = IniReadReal ('topt', speciesDB, section = ToString(i)) species (i) % theta_fswc = IniReadReal ('phi-theta', speciesDB, section = ToString(i)) species (i) % theta_fvpd = IniReadReal ('phi-ea', speciesDB, section = ToString(i)) species (i) % fpra = IniReadReal ('fpra', speciesDB, section = ToString(i)) species (i) % fprn = IniReadReal ('fprn', speciesDB, section = ToString(i)) species (i) % spra = IniReadReal ('spra', speciesDB, section = ToString(i)) species (i) % sprn = IniReadReal ('sprn', speciesDB, section = ToString(i)) species (i) % ltr = IniReadReal ('leaf-turnover', speciesDB, section = ToString(i)) species (i) % rtr = IniReadReal ('root-turnover', speciesDB, section = ToString(i)) species (i) % tcold_leaf = IniReadReal ('tcold-leaf', speciesDB, section = ToString(i)) species (i) % sla = IniReadReal ('sla', speciesDB, section = ToString(i)) species (i) % hdmin = IniReadReal ('hdmin', speciesDB, section = ToString(i)) species (i) % hdmax = IniReadReal ('hdmax', speciesDB, section = ToString(i)) species (i) % albedo = IniReadReal ('albedo', speciesDB, section = ToString(i)) species (i) % laimax = IniReadReal ('laimax', speciesDB, section = ToString(i)) species (i) % canopymax = IniReadReal ('canopymax', speciesDB, section = ToString(i)) species (i) % wood_density = IniReadReal ('wood-density', speciesDB, section = ToString(i)) species (i) % ms = IniReadReal ('ms', speciesDB, section = ToString(i)) species (i) % mr = IniReadReal ('mr', speciesDB, section = ToString(i)) species (i) % mf = IniReadReal ('mf', speciesDB, section = ToString(i)) species (i) % wSx1000 = IniReadReal ('wSx1000', speciesDB, section = ToString(i)) END DO ! species DB is complete. Deallocate ini database CALL IniClose (speciesDB) RETURN END SUBROUTINE ReadSpecies !============================================================================== !| Description: ! set initial condition SUBROUTINE SetIC & ! (inifile) IMPLICIT NONE CHARACTER (LEN = *), INTENT(in) :: inifile !local declarations: TYPE (IniList) :: icDB INTEGER (KIND = short) :: ncohorts TYPE (grid_integer) :: cohort_map INTEGER (KIND = short) :: i, j, k, l INTEGER (KIND = short) :: ispecies REAL (KIND = float) :: dbh TYPE (PlantsCohort), POINTER :: cohort !------------------------------end of declarations---------------------------- ! open and read configuration file CALL IniOpen (inifile, icDB) !number of cohorts ncohorts = IniReadInt ('cohorts', icDB) !starting year year_new = IniReadInt ('year', icDB) year_prev = IniReadInt ('year', icDB) !cohort map CALL GridByIni (icDB, cohort_map, section = 'cohorts-map') DO k = 1, count_stands !aggiungere configurazione ricorsiva quando ci sono più coorti per cella !al momento solo una coorte per stand (cella) cohort => forest (k) % first DO l = 1, forest (k) % lenght ispecies = IniReadInt ('species', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % species = species (ispecies) cohort % age = IniReadInt ('age', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % density = IniReadReal ('density', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % height = IniReadReal ('height', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % dbh = IniReadReal ('dbh', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % lai = IniReadReal ('lai', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % mass_stem = IniReadReal ('stem-biomass', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % mass_root = IniReadReal ('root-biomass', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % mass_leaf = IniReadReal ('leaf-biomass', icDB, section = ToString (cohort_map % mat (forest (k) % i, forest (k) % j ) ) ) cohort % mass_total = cohort % mass_stem + cohort % mass_root + cohort % mass_leaf cohort % crown_diameter = CrownDiameter ( dbh = cohort % dbh, den = cohort % density, & denmin = species (ispecies) % denmin, & denmax = species (ispecies) % denmax, & dbhdcmin = species (ispecies) % dbhdcmin, & dbhdcmax = species (ispecies) % dbhdcmax ) cohort % canopy_cover = CanopyCover ( cohort % crown_diameter, cohort % density ) cohort => cohort % next END DO END DO ! Deallocate memory CALL IniClose (icDB) RETURN END SUBROUTINE SetIC !============================================================================== !| Description: ! Computes the absorbed photosynthetically active radiation. ! The amount of light available for the plant that is going to define ! its growth rate. The light that the plant could absorb is determined according ! to the canopy cover. The amount of absorbed photosynthetically active radiation ! is usually computed following Lambert-Beer law . ! ! References: ! ! McCree, K.J.: Test of current definitions of photosynthetically active ! radiation against leaf photosynthesis data. Agr. Forest. Meteorol, ! 10 , 443-453, 1972. FUNCTION AparCalc & ! (rad, lai, k, alb, dt) & ! RESULT (apar) IMPLICIT NONE !arguments with intent(in): REAL (KIND = float), INTENT(IN) :: rad !! shortwave direct radiation (w/m2) REAL (KIND = float), INTENT(IN) :: lai !! leaf area index (m2/m2) REAL (KIND = float), INTENT(IN) :: k !! extinction coefficient for absorption of PAR by canopy REAL (KIND = float), INTENT(IN) :: alb !! plant albedo (0-1) INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s) !local declarations: REAL (KIND = float) :: apar REAL (KIND = float) :: par0 !!photosynthetically active radiation (molPAR m-2) REAL (KIND = float), PARAMETER :: radToPAR = 4.57 !conversion factor 1 W m-2 = 4.57 µmol m-2 s-1 (McCree, 1972) !------------------------------end of declarations---------------------------- !devo considerare l'albedo cioè è una short wave net radiation ? o quella incidente ? par0 = rad * ( 1. - alb ) * radToPAR * dt ! (µmol m-2) apar = par0 * ( 1. - EXP ( -k * lai) ) / 1000000. ! (molPAR m-2) RETURN END FUNCTION AparCalc !============================================================================== !| Description: ! update leaf biomass and leaf area index ! ! References: ! !Maneta, M. P., and N. L. Silverman, 2013: A spatially distributed model to ! simulate water, energy, and vegetation dynamics using information from ! regional climate models. Earth Interact., 17 ! doi:10.1175/2012EI000472.1. ! ! Arora, V. K., and G. J. Boer, 2005: A parameterization of leaf phenology ! for the terrestrial ecosystem component of climate models. ! Global Change Biol., 11, 39–59. ! SUBROUTINE GrowLeaf & ! (npp, af, Ta, Tcold, swc, swp, sfc, tr, sla, mleaf, lai) IMPLICIT NONE !arguments with intent(in): REAL (KIND = float), INTENT(IN) :: npp !!net primary production [t/ha] REAL (KIND = float), INTENT(IN) :: af !!allocation factor [0-1] REAL (KIND = float), INTENT(IN) :: Ta !! air temperature [°C] REAL (KIND = float), INTENT(IN) :: Tcold !! temperature threshold that accelerates leaf turnover (°C) REAL (KIND = float), INTENT(IN) :: swc !! soil water content [m3/m3] REAL (KIND = float), INTENT(IN) :: swp !! soil wilting point [m3/m3] REAL (KIND = float), INTENT(IN) :: sfc !! soil field capacity [m3/m3] REAL (KIND = float), INTENT(IN) :: tr !! leaf turnover rate [s-1] REAL (KIND = float), INTENT(IN) :: sla !! specific leaf area [m2/Kg] !arguments with intent(inout): REAL (KIND = float), INTENT(INOUT) :: mleaf !!mass of leaf [t/ha] REAL (KIND = float), INTENT(INOUT) :: lai !! leaf area index [m2/m2] !local declarations: REAL (KIND = float) :: deltaMleaf !! leaf mass increment REAL (KIND = float) :: swc_stress !! soil water content stress that affects leaf increment [0-1] REAL (KIND = float) :: temp_stress !! temperature stress that affects leaf increment [0-1] REAL (KIND = float) :: actual_tr !! actual leaf tirnover rate with temperature and soil water content stresses !------------------------------end of declarations---------------------------- !compute temperature stress IF (Ta >= Tcold) THEN temp_stress = 1. ELSE IF ((Tcold - 5.) < Ta < Tcold) THEN temp_stress = (Ta - Tcold - 5.) / 5. ELSE IF (Ta <= (Tcold - 5.)) THEN temp_stress = 0. END IF !temp_stress = MAX (0., MIN (1., ((Ta - (Tcold - 5.) / 5. ) ) )) !compute hydrologic stress swc_stress = MAX (0., MIN (1., (swc - swp) / (sfc - swp) ) ) !compute actual leaf turnover rate: base turnover + 0.001 hydrologic stress + temperature stress actual_tr = tr * ( 1. + 0.001 * (1. - swc_stress)**3. + (1. - temp_stress)**3. ) !leaf mass increment deltaMleaf = af * npp !update leaf biomass mleaf = mleaf + deltaMleaf !update leaf area index lai = lai + (deltaMleaf * sla/ 10.) - (actual_tr * dtPlants) *LAI !10 = conversion factor 1000 [kg/t] / 10000 [m2/hectare] !check lower boundary IF (mleaf < 0.) THEN mleaf = 0. END IF IF (lai < 0.) THEN lai = 0. END IF RETURN END SUBROUTINE GrowLeaf !============================================================================== !| Description: ! update DBH and height tree according to ech2o model by Maneta ! SUBROUTINE GrowDBHech2o & ! (cc, hdmin, hdmax, ws, dws, DBH, height, tree_density, wood_density, age, maxage) IMPLICIT NONE !arguments with intent(in): REAL (KIND = float), INTENT(IN) :: cc !! canopy cover REAL (KIND = float), INTENT(IN) :: hdmin !! H/D ratio in carbon partitioning for low density REAL (KIND = float), INTENT(IN) :: hdmax !! H/D ratio in carbon partitioning for high density REAL (KIND = float), INTENT(IN) :: ws !! stem biomass (t/ha) REAL (KIND = float), INTENT(IN) :: dws !! stem biomass increment (t/ha) REAL (KIND = float), INTENT(IN) :: tree_density !! tree density (trees/ha) REAL (KIND = float), INTENT(IN) :: wood_density !! wood density (kg/m3) REAL (KIND = float), INTENT(IN) :: age !! tree age (year) REAL (KIND = float), INTENT(IN) :: maxage !!tree max age (year) !arguments with intent inout REAL (KIND = float), INTENT(INOUT) :: DBH !diameter at brest height (cm) REAL (KIND = float), INTENT(INOUT) :: height !tree height (m) !local declarations REAL (KIND = float) :: delta_tree_mass !!tree mass increment (kg/tree) REAL (KIND = float) :: fhd !! grow factor that depends on the height-to-diameter ratio REAL (KIND = float) :: dDBH !DBH increment (cm) REAL (KIND = float) :: dheight !height increment (m) !-------------------------------end of declarations---------------------------- !compute mass increment per tree delta_tree_mass = dws / tree_density * 1000. ! mass increment per tree (kg) !set grow factor IF (height / DBH >= hdmin .AND. cc < 0.95) THEN fhd = hdmin ELSE IF (height / DBH <= hdmax .AND. cc >= 0.95 .AND. age < maxage / 2) THEN fhd = hdmin ELSE IF (height / DBH <= hdmax .AND. cc >= 0.95 .AND. age > maxage / 2 ) THEN fhd = hdmax ELSE IF (height / DBH < hdmin) THEN fhd = hdmax ELSE IF (height / DBH > hdmax) THEN fhd = 0.5 * hdmin ELSE IF ( age > 0.7 * maxage) THEN fhd = 0. END IF !compute DBH increment dDBH = 4. * delta_tree_mass / ( wood_density * pi * (DBH / 100.)**2. * ( 2. * height / (DBH / 100.) + fhd * 100. ) ) * 100. !compute tree height increment dheight = dDBH * fhd !update dbh and height tree DBH = DBH + dDBH height = height + dheight !write(*,*) dDBH, dheight !pause RETURN END SUBROUTINE GrowDBHech2o !============================================================================== !| Description: ! update DBH and height tree every new year ! SUBROUTINE GrowDBH & ! (cc, hdmin, hdmax, ws, dws, DBH, height) IMPLICIT NONE !arguments with intent(in): REAL (KIND = float), INTENT(IN) :: cc !! canopy cover REAL (KIND = float), INTENT(IN) :: hdmin !! H/D ratio in carbon partitioning for low density REAL (KIND = float), INTENT(IN) :: hdmax !! H/D ratio in carbon partitioning for high density REAL (KIND = float), INTENT(IN) :: ws !! stem biomass (t/ha) REAL (KIND = float), INTENT(IN) :: dws !! stem biomass increment (t/ha) !arguments with intent inout REAL (KIND = float), INTENT(INOUT) :: DBH !diameter at brest height (cm) REAL (KIND = float), INTENT(INOUT) :: height !tree height (m) !local declarations REAL (KIND = float) :: hdeff REAL (KIND = float) :: dDBH !DBH increment (cm) REAL (KIND = float) :: dheight !height increment (m) !-----------------------------end of declarations------------------------------ IF ( cc <= 0.95) THEN !low density hdeff = hdmin + (hdmax - hdmin) * cc ELSE ! high density hdeff = hdmax END IF !compute height increment !dheight = hdeff / ( hdeff / height + 200. / dbh ) * dws / ws dheight = hdeff / ( hdeff / height) + (200. / dbh ) * dws / ws !compute DBH increment dDBH = ( height / hdeff ) + ( 200. / dbh ) * ( dws / ws ) !update dbh and height tree DBH = DBH + dDBH height = height + dheight !write(*,*) dheight, hdeff, dws, ws !pause RETURN END SUBROUTINE GrowDBH !============================================================================== !| Description: ! fill in a grid with state variable values in forest stands ! SUBROUTINE ForestToGrid & ! (grid, statevar) IMPLICIT NONE !arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: statevar !arguments with intent inout TYPE (grid_real), INTENT (INOUT) :: grid !local declarations INTEGER (KIND = short) :: i, j, k !-----------------------------end of declarations------------------------------ DO k = 1, count_stands i = forest (k) % i j = forest (k) % j SELECT CASE (statevar) CASE ('lai') !leaf area index grid % mat (i,j) = forest (k) % first % lai CASE ('fv') !canopy cover grid % mat (i,j) = forest (k) % first % canopy_cover CASE ('gpp') !GPP grid % mat (i,j) = forest (k) % first % gpp * & CellArea (grid,i,j) / hectare CASE ('npp') !NPP grid % mat (i,j) = forest (k) % first % npp * & CellArea (grid,i,j) / hectare CASE ('root') !root mass grid % mat (i,j) = forest (k) % first % mass_root * & CellArea (grid,i,j) / hectare CASE ('stem') !stem mass grid % mat (i,j) = forest (k) % first % mass_stem * & CellArea (grid,i,j) / hectare CASE ('leaf') !foliage mass grid % mat (i,j) = forest (k) % first % mass_leaf * & CellArea (grid,i,j) / hectare CASE ('stemyield') !stem mass grid % mat (i,j) = forest (k) % first % stem_yield * & CellArea (grid,i,j) / hectare CASE ('dbh') !diameter at brest heigth (cm) grid % mat (i,j) = forest (k) % first % dbh CASE ('height') !tree height (m) grid % mat (i,j) = forest (k) % first % height CASE ('density') !tree density (trees/hectare) grid % mat (i,j) = forest (k) % first % density END SELECT END DO RETURN END SUBROUTINE ForestToGrid END MODULE Plants